function varargout=compute_MBC_shock4(varargin)
% This function computes MBC shock (Main Business Cycle) denoted as q
% input
% Compared to compute_MBC_shock, this allows multiple variables 
% varargin{1}=out: VAR result
% varargin{2}=ps: parameters
% varargin{3}=opt: option
% 
% output is rs including q
% 
% 2020/6/29 copied from compute_MBC_shock3
% use bins between [0,2*pi]

%--------------------------------------------------------------------------
% input
out=varargin{1};
ps=varargin{2};
opt=varargin{3};
%--------------------------------------------------------------------------
% compute H 
H=zeros(ps.N,ps.N); % integral of spectral density D(w) associated with q
D_vec=zeros(ps.ngrid,ps.N,ps.N); % D(w)
D_all_vec=zeros(ps.ngrid,1); % spectral density at each w

qDq_vec=zeros(ps.ngrid,ps.N,ps.N);

bin_vec=[];%zeros(ps.nbin,1); % bin values, not used anymore


for i=1:ps.ngrid % grids [0,2*pi]

% compute D(w)    
tmp=compute_Dw4(ps.grid(i),out,ps,opt);

% to compute q'D(w)q (spectral density associated with q at w)
qDq_vec(i,:,:)=tmp.D; % D(w) at w

% q shock, matrix D(w) in [w1,w2] with width
D_vec(i,:,:)=tmp.D*ps.igrid(i)*(2*pi/ps.ngrid); % D(w)*width,0 if grid outside [w1,w2]

% all shocks (scalar, spectral density at w)
D_all_vec(i)=tmp.D_all;

H=H+squeeze(D_vec(i,:,:)); % integral of D(w) in [w1,w2]

end


% compute eigenvector

% only largest eigenvalue, faster
% [eig_vec,eig_val]=eigs(H,1); % 1st largest eigenvalue and eigenvector q 
% q=eig_vec(:,1); % 1st largest eigenvector

% if want to get all eigenvalues, slower
[eig_vec,eig_val]=eig(H); % all eigenvalues
[d,ind]=sort(diag((eig_val)),'descend'); % largest to smallest order
eig_val=eig_val(ind,ind); % reorder eigenvalue
eig_vec=eig_vec(:,ind); % reorder eigenvector
q=eig_vec(:,1); % 1st largest eigenvalue 


% normalize the eigenvector
q=q/norm(q); % |q|=1 

% if q(ps.k)<0 % q(k)>0 : k-th row of q is positive
%     q=-q;
% end
S_tilda=chol(out.SIGMA,'lower');
tmp=S_tilda*q;

if tmp(ps.k(1))<0
    q=-q;
end
% compute D(w) for q
qDq=zeros(ps.ngrid,1);
for i=1:ps.ngrid
    qDq(i)=q'*squeeze(qDq_vec(i,:,:))*q;
end

%--------------------------------------------------------------------------
% output
rs.q=q;
rs.eig_vec=eig_vec;
rs.eig_val=eig_val;
rs.D_all_vec=D_all_vec;
rs.bin_vec=bin_vec;
rs.qDq=qDq;
varargout{1}=rs;




